#!/bin/bash -e
function info() {
echo Usage: `basename $0` '[-st] r1 r2'
exit 2
}

while getopts  ":p:f:s:t:" opt; do
	case  $opt  in
		p) out_prefix=$OPTARG;;
		f) suffix=$OPTARG;;
        s) sample_name=$OPTARG;;
        t) run_continue_dir=$OPTARG;;
		*) info;;
	esac
done
shift $(($OPTIND - 1))

if [ $# -lt 2 ]; then info; fi



# makedir=${makedir:=T}

. $var

mnt=mnt.sh
path_abs=path_abs.sh
makedir=makedir.sh
sam2b=sam2b_samtools.sh
merge_bam=merge_bam.sh
summ=summ.sh
summ_gatk=summ_gatk.sh

samtools_path=/mnt/ilustre/app/rna/bin
java=$tools_path/jre1.6.0_45/bin/java
java_run1="$java $java_mem -jar $tools_path/jars"
# export CLASSPATH=.:$JAVA_HOME/lib/dt.jar:$JAVA_HOME/lib/tools.jar


r1=$(path_abs.sh $1)
r2=$(path_abs.sh $2)


if test -z "$run_continue_dir"; then
run_continue_dir=`$makedir -p$sample_name` && cd $run_continue_dir
else 
cd $run_continue_dir
fi

$mnt $$ `basename $0 .sh`.`basename $run_continue_dir` & # `basename $run_continue_dir` for strip / if exist in dir's end.


echo;echo $r1 $r2
adapters_file=$tools_path/Trimmomatic-0.36/adapters/TruSeq3-PE.fa
primary_bed=$data_path/intervals/1/primary.bed 
capture_bed=$data_path/intervals/1/capture.bed
ref_genome_size=ref_genome.size.txt


# Determine methylation percentage using BSMAP
# support samtools 0.19.
# very memory eatting, min 30g
# echo;echo Determine methylation percentage 
# methratio.py  -d $ref_genome \
# -s $samtools_path \
# -m 1 -z -i skip \
# -o $out_prefix.methylation_results.txt \
# $out_prefix.clipped.bam 

# echo;echo Determine methylation percentage 
# methratio.py  -d $ref_genome \
# -s samtools -m 1 -z -i skip \
# -c $sample_name \
# -o $out_prefix.$sample_name.methylation_results.txt \
# $out_prefix.clipped.bam 


# Determine bisulfite conversion efficiency using BSMAP
echo;echo Base Quality Recalibration first step
$java_run1/bissnp -R $ref_genome \
-I $out_prefix.clipped.bam \
-T BisulfiteCountCovariates \
-cov ReadGroupCovariate \
-cov QualityScoreCovariate \
-cov CycleCovariate \
-recalFile $out_prefix.recalFile_before.csv \
-knownSites $data_path/ncbi/dbsnp/All_20150605.vcf
# -nt 4


echo;echo Base Quality Recalibration second step
$java_run1/bissnp -R $ref_genome \
-I $out_prefix.clipped.bam \
-o $out_prefix.recal.bam \
-T BisulfiteTableRecalibration \
-recalFile $out_prefix.recalFile_before.csv \
-maxQ 40


echo;echo Combined SNP/methylation calling 
$java_run1/bissnp -R $ref_genome \
-I $out_prefix.recal.bam \
-T BisulfiteGenotyper \
-vfn1 $out_prefix.cpg.raw.vcf  \
-vfn2 $out_prefix.snp.raw.vcf  \
-L $capture_bed \
-stand_call_conf 20 \
-stand_emit_conf 0  \
-mmq 30 \
-mbq 0  \
-nt 4
# -D genome.snps.vcf \


echo;echo Sort VCF files 
sortByRefAndCor.pl --k 1 --c 2 \
$out_prefix.snp.raw.vcf \
$ref_genome.fai \
> $out_prefix.snp.raw.sorted.vcf 

sortByRefAndCor.pl --k 1 --c 2 \
$out_prefix.cpg.raw.vcf \
$ref_genome.fai \
> $out_prefix.cpg.raw.sorted.vcf 


 
echo;echo Filter SNP/methylation calls  
$java_run1/bissnp -R $ref_genome \
-T VCFpostprocess \
-oldVcf $out_prefix.snp.raw.sorted.vcf \
-newVcf $out_prefix.snp.filtered.vcf  \
-snpVcf $out_prefix.snp.raw.sorted.vcf \
-o $out_prefix.snp.filter.summary.txt  
 
$java_run1/bissnp -R $ref_genome \
-T VCFpostprocess \
-oldVcf $out_prefix.cpg.raw.sorted.vcf \
-newVcf $out_prefix.cpg.filtered.vcf  \
-snpVcf $out_prefix.snp.raw.sorted.vcf \
-o $out_prefix.cpg.filter.summary.txt

echo;echo Convert VCF to BED file
vcf2bed6plus2.pl  $out_prefix.snp.filtered.vcf
vcf2bed6plus2.pl  $out_prefix.cpg.filtered.vcf

. $cmd_done